principals  and  officers 


. , . computers  to  the 

had\/PY  President  ■ e application  . . ;n  the  areas 

COY  P-  HARvt  , experience  w and  particularly  software 

c——-. 

* crane  senior  Vice  President  ^ ^ sysKrns^nS 

M'CHA0nOon.  .n«r-.ro-. 

University.  He*  , sc-,ence  from  P rams  and  his  text  • ^qcj  extensive 

*•  ^SS'lSod  to*  ““"iSteSndard  » *' airline  Indes....* 

**  p°pe" 


in  w w ^ 

DONALD  L.  ^3^^^ 

gSSsBS-  ““ 


*•  naci  ies in  industrial  ana 
jlting  capacities 

_ .»nPFY  Principal  , e areas  of  'ieJa't^ ' ^ allocation. 

^ has  had  extensive  exP*^enryB  management,  °nd « president,  Stan- 

<«0  K‘“"h  '"*  ’ ..errors  o«  *» 


£&“='- or.  P-!-  — - 

ANALYSIS  CORPOPATlOPr*  Independenr  auditors 
CONTROL  ANALYSIS 
Company- 


_ MCUtMNMV 

lilt  *»«•  tKltu  X 

rs#  toft  toctw  □ 

ttimtoircn  □ 

jtnintaxa 

«v  Q.Yv  ^-Ac 

t -I  swa 


3$;  ■ 


3 2*  f ' 


Eb 


/o 


CONTROL  ANALYSIS  CORPORATION 
800  Welch  Road 
Palo  Alto,  California 


/Technical  j/leport^No.  86-24 
/TH  November  1-977  f 


ik^ 1 


CONFIDENCE  .LIMITS 
FOR  THE  ^OPTIMUM 


SIMULATION 
RESPONSE  FUNCTION  s 


lOl 


chael  A./ Crane  and  David  W.jRobinson  / 


D D C 

NO V 27  1978 


nl  ftMt- work  was  .supported  by  the  Office  of  Naval  Research  under  contract 
/O >10001 4-  72-C-0O86/ NR-047- 106). 

Reproduction  in  whole  or  in  part  is  permitted  for  any  purpose  of  the 
United  States  Government. 

j^^UTjQN  STATEMENT  A ’ | 

Approved  fa*  public  2llMle;  i 
— Piittibufton  Unlimited  I 


TABLE  OF  CONTENTS 


Paae 

1.  INTRODUCTION  AND  SUMMARY  1 

2.  BASIC  RESULTS  4 

3.  NUMERICAL  EXAMPLE 11 

4.  SOME  EXTENSIONS I5 

4. 1 Additional  Observations  15 

4.2  Interval  Width  Selection  16 

4.3  Non- Concave  Functions  20 

5.  REFERENCES 26 


i 


1. 


INTRODUCTION  AND  SUMMARY 


oThis  paper  deals  with  stochastic  simulation  in  which  the  system  being 
studied  can  be  controlled  to  some  extent  by  setting  the  value  of  an  input 
parameter.  The  system  response  for  a given  input  will  vary  from  simulation 
run  to  simulation  run  due  to  uncontrolled  and  random  factors;  we  thus  con- 
sider the  average  or  expected  response.  This  is  usually  called  the  simula- 
tion response  surface.  Since  our  work  here  involves  only  a single  variable 
the  term  response  curve  might  be  more  appropriate. ^ 

1 — The  simulation  "response"  will  generally  be  identified  with  some  measure 
of  system  effectiveness  and  the  analyst  will  want  to  maximize  (or  possibly 
minimize)  the  response  by  choosing  the  optimum  value  for  the  input  parameter. | 
More  formally,  if  we  denote  the  input  parameter  by  X and  the  expected 
response  by  g(X)  , we  wish  to  solve  the  problem 

maximize  g(x)  ( 1) 

The  main  difficulty  in  solving  ( l)  is  that  g(x)  can  be  evaluated  only 
approximately  due  to  the  random  factors  in  the  simulation.  If  a simulation 
is  carried  out  with  the  input  parameter  X fixed  at,  say,  a^  , the  output 
from  the  run  will  be  a confidence  interval  for  which 


PrfAj1  < g(a]L)  < A~  ) = 1 - 71  , 


1 2 

where  y^  is  a user-specified  value  and  [A^  , A“  ] is  the  confidence 
interval.  Such  intervals  depend  on  some  assumed  probability  distribution 
(often  normal)  for  the  simulation  outputs.  We  do  not  consider  this  under- 
lying distribution  more  explicitly  in  this  paper;  the  reader  may  refer  to 
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L 


i 


[2]  or  [3]  for  some  examples.  We  do  require  that  intervals  of  the  form  (2) 
be  found  for  each  simulation  run  that  is  carried  out. 

The  basic  result  of  this  paper  requires  as  given  a set  of  confidence 
intervals  (A^)  for  the  response  surface  values  at  fixed  input  parameter 
values  la^)  . (There  are  a variety  of  approaches  for  obtaining  these  con- 
fidence intervals,  e.g.  batch  means,  independent  replications,  autoregressive 

schemes,  regenerative  approach,  etc.)  These  are  used  to  obtain  confidence 

* * 

intervals  for  the  solution  X to  (1),  the  optimum  value  g(X  ) and  a 

* * 

joint  confidence  region  for  the  point  (X  , g(X  ))  , in  the  case  where  the 
response  surface  is  quadratic,  i.e. 


g(x)  = zl  + z2X  + z^> 


The  technique  applied  here  is  related  to  a confidence  band  methodology 

developed  in  a number  of  previous  papers  ([2]  - [4]).  The  basic  idea  is  to 

find  the  locus  of  the  optima  of  all  possible  quadratic  functions  which  pass 

simultaneously  through  each  of  the  confidence  intervals.  The  result  (under 

some  mild  convexity  assumptions)  is  a compact  region  whose  x-axis  dimensions 

•* 

give  the  required  confidence  interval  for  X . The  probability  coverage 
of  the  final  region  is  found  by  multiplying  together  the  probabilities  of 
the  individual  response  surface  confidence  intervals. 

The  main  drawback  to  our  method  is  that  the  probability  coverage  of 
the  confidence  interval  for  the  optimum  may  be  quite  low  when  many  different 
parameter  settings  are  used.  By  restricting  ourselves  to  a quadratic  re- 
sponse surface  with  a single  input  parameter,  then,  we  may  obtain  results 


with  but  three  or  four  observations  while  keeping  the  probability  coverage 


high.  Moreover,  the  experience  of  workers  in  nonlinear  programming  and 
approximation  theory  indicates  that  a quadratic  approximation  is  in  many 
cases  an  entirely  adequate  representation  for  a smooth  function  over  a 
restricted  range.  Thus,  the  assumed  form  (3)  for  g(x)  is  not  necessarily 
restrictive. 

The  major  strength  of  our  technique  is  that  it  does  not  depend  on  the 
distribution  of  the  simulation  outputs  in  any  way.  In  fact,  probabilities 
are  used  only  to  find  the  probability  of  the  final  region  by  a simple  multi- 
plication. There  is  no  need  for  regression,  for  equal  variance  assumptions 
or  for  distribution  theory  in  analyzing  simulation  outputs. 

Most  other  approaches  to  finding  the  optimum  of  a simulation  response 
function  involve  carrying  out  a search  over  the  parameter  space;  see  [5] 
for  a review  and  evaluation  of  work  in  this  area.  In  our  method  no  explicit 
attempt  is  made  to  determine  the  parameter  setting  for  the  next  simulation 
run;  rather,  the  emphasis  is  on  using  the  information  already  obtained  to 
determine  where  the  optimum  might  be.  The  best  use  of  the  present  work, 
then, would  be  to  carry  out  an  analysis  of  the  final  results  of  a simulation, 
perhaps  using  only  the  last  3”5  parameter  settings. 

The  remainder  of  the  paper  is  organized  as  follows:  Section  2 presents 
the  basic  results  and  Section  3 consists  of  a numerical  example.  The  final 
section  reports  on  some  straightforward  extensions  to  the  basic  results  and 
applies  them  to  the  numerical  example. 


2. 


BASIC  RESULTS 


Suppose  that  the  response  surface  g(  X)  is  as  given  by  (3)  and  that 
simulation  runs  have  been  performed  for  the  input  parameter  settings 
X = a^  > 1 < i < n y where  a^  < . . . < a^  , n > 3 . Suppose 

further  that  the  simulation  results  are 

1 2 

Ai  - 8^ai^  - Ai  with  probability  1 - 7^  , 

A^  < g(a2)  < A|  with  probability  1 - 7g  , (4) 

• • • 

1 2 

A < g(  a ) < A with  probability  1 - 7 

n — n'  — n J /n 

Assuming  the  simulation  runs  are  independent,  g(x)  will  satisfy  all  of  the 
inequalities  (4)  with  probability  (l  - 7 ^) ( 1 - 7 ) * ' * (1  - 7^)  = 1-7 
(Alternatively,  one  might  deal  with  simulations  with  dependent  results  as 
long  as  the  probability  1 - 7 may  be  found;  this  seems  unlikely  in  practice, 
however,  so  we  will  not  consider  it  further  here.) 

With  probability  1 - 7 , then,  g (x)  must  pass  through  the  points 

(a^  , y^)  , i = 1,  . . . n , where  A,,  < y^  < A(~  . Consider  the  vector 
y = (yv  yo,  . . . , yn)  ; if  three  of  the  components  of  y are  fixed,  the 
corresponding  value  of  z - ( z^,  z,„  z .)  from  (3)  may  be  found  and  the 
remaining  components  of  y evaluated.  Assuming  without  loss  of  generality 
that  the  three  fixed  components  are  y^,  y ,,  y . we  have 
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z = 


(yl  * y2  ' y5) 


(a5“  ai)(a2“  ai)(a3-  a2) 


(V  a2)a2a3  a2"Sj 


(a2“ai)aia£ 


Va2 


(al"a5)ala3  Val  al"a3 


2 2 

a,  - a - a 

12  2 l 


(5) 


Thus  (5)  gives  the  coefficients  of  a quadratic  function  passing  through  the 
points  (a^,  y£)  > i = 1,  2,  3 • Note  that  the  components  of  z in  (5)  are 
simple  linear  functions  of  the  components  of  y . We  may  thus  refer  to 
MZj(y)"  , say,  to  indicate  the  result  of  (5)  for  a specific  choice  of  the 
y vector. 

When  there  are  more  than  three  simulation  runs  it  is  likely  that  not  all 
possible  choices  for  the  components  of  y will  be  admissible  since  there  may 
be  no  quadratic  function  passing  through  the  specified  points.  We  thus  require 
that  y be  a member  of  the  constraint  set  Y , where 


= (y f Aj1  < yi  < A.2  ; i = 1,  2,  3 ; 


(6) 


i < zL(y)  + z )(y)ai  + z (y)a^  < ; 3 < i < n) 


Note  that  all  of  the  constraints  in  Y are  linear. 

The  optimum  of  the  quadratic  function  (3)  occurs  at  the  parameter  setting 


X = - z2/2z5  . (7) 

This  will  be  a maximum  as  long  as  z < 0 , i.e.  the  response  surface  is 

strictly  concave.  We  assume  for  the  present  that,  in  fact,  z„  < 0 ; dealing 
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with  linear  or  convex  response  surfaces  is  discussed  more  fully  in  Section 


The  X 
z and  z_ 

2 3 


value  given  by  (7)  may  be  viewed  as  a function  of  y just  as 
are.  Substituting  from  (5)  and  simplifying  yields 


x(y) 


/ 2 2\  / 2 2»  ,2  2S 

- ag  )y1-  (a^  - ax  )y£  +(ag  - &1  )y^ 

2(a3-  a2)yx-  2(a5-  a^y^  2(ag-  a^ 


(8) 


Thus  (8)  gives  the  optimum  X of  a quadratic  function  passing  through  the 
points  (a^,  y^)  , i = 1,  2,  3 • Since  we  have  assumed  z^  < 0 , X(y) 
is  a continuous  bounded  function  over  the  closed  compact  set  Y . Thus, 

A 

as  a function  of  yeY  , X(y)  takes  on  a maximum  (x)  and  a minimum 
(x)  and  all  intervening  values. 

Now,  given  that  g(x)  is  in  fact  quadratic  and  that  it  satisfies  (4) 
(i.e.,  is  contained  by  all  of  the  simulation  output  confidence  intervals) 
the  optimum  input  parameter  must  lie  in  the  closed  interval  [X,  X]  , 

— /\  /S 

where  X = sup  X(y)  and  X = inf  X(y)  . Since  the  probability  that  (4) 
ycY  ycY 

holds  is  just  I-7  , we  have 


Pr{x  < X*  < X)  < I-7 


(9) 


This  is  the  desired  confidence  interval  for  X . (The  inequality  results 
from  the  observation  that  the  optimal  parameter  may  lie  in  the  given  interval 
even  in  cases  where  the  quadratic  does  not  lie  within  the  simulation  output 
confidence  interval.) 
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Finding  X and  ^ is  quite  straightforward.  Examining  the  gradient 
vector  found  by  differentiating  (0): 


VyX(y) 


2f(a^- 


(a2-a1)(a,-a;J(ay-a1) 

a2)yr(a3'ai)y2+(a2-ai)y3r" 


y3 


shows  that  the  fraction  will  always  be  positive,  since  the  denominator  is 

2 

proportional  to  zv  . Thus,  when  two  of  the  first  three  components  of  y 

A 

are  fixed  the  direction  of  change  of  X(y)  is  constant.  (The  rate  of 

A 

change  varies  but  not  the  direction.)  X(y)  may  thus  be  increased,  say,  by 
changing  the  unfixed  y component  until  one  of  the  constraints  in  (6) 
becomes  binding;  this  shows  that  both  _X  and  X occur  at  extreme  points 
of  the  set  Y , i.e.  three  of  the  simulation  output  confidence  interval 
inequalities  (4)  will  be  binding. 

Now  suppose  that  X is  fixed  at  some  point  f_X,  X]  . We  know  that 
there  is  at  least  one  y vector  such  that  X( y)  = X , i.e.  the  quadratic 
function  passing  through  ( a^  , y^)  , 1 < i < n , has  its  optimum  at  X 

In  general  there  will  be  a set  of  such  vectors  and  the  maximum  response 
value  g(x)  will  vary  as  different  y values  are  selected.  Denote  by 
g( X ; y)  this  maximum  and  define  the  functions 


_s(  > ) = inf  g(  X ; y)  ( 10) 

X(y)  = X 

y ■ Y 
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r 


g(  A ) = sup  g( ; ; y) 

x(y)  = A 

y c Y 

These  functions  then  depict  a joint  confidence  region  for  the  point 
(A*,  g(X*))  , i.e. 

Pr((X*,  g(x*))  r G)  > i - 7 


where 


G = C(X  , g)  | X < X < X ; g(x)  < g<  g(x) ) . 

This  follows  just  as  for  (9). 

Now  (10)  and  (11)  are  linear  programs  for  fixed  X . We  show  this 
as  follows:  the  first  constraint  (x(y)  = X)  may  be  rewritten  as 
-z0(y)  = 2Xz^( y)  according  to  ('/);  simplifying  using  (5)  yields 


(ah“  ap^ah  +ap"  2X)  (a0-  +a1  - 2X) 

y = ^ ^ — y + — — = y 

(a^-  a1)(a5  +aL-  2k)  ( - a1)(a5  +aL-  2X)  5 


(n) 


(12) 


(13) 


(13)  shows  that  as  y^  and  y.  change,  y(j  shifts  so  as  to  keep  the  opti 

fixed  at  X The  value  of  g(X  ; y)  is  given  by  z.(y)  - z (y)X‘~  ; 

1 3 

substituting  from  {'})  and  ( I3)  gives 


imum 


l 


g(x  ; y)  = 


(y*)‘ 


y i - 

(y  ax)(a^  +ai  - 2X)  (a  - a^a^  2X) 


(x-a^ 


(14) 


The  solution  to  ( 10 ) and  (11)  will  then  have  two  components  of  y (say,  y^ 
and  y ) at  their  bounds,  with  the  third  component  y^  determined  by  the 

A 

optimality  constant  X(y)  = X , or  (lj). 

This  analysis  must  be  repeated  when  a^  + a^  - 2X  = 0 , i.e.  X is 

at  the  midpoint  of  [a,  , a ] . (Note  that  this  point  may  not  be  in 

1 J 

[X,  X].)  In  this  case,  we  obtain 


yl  = 


(15) 


There  are  two  alternatives  for  computing  the  bound  functions.  If  this 
is  to  be  done  by  "brute  force"  (either  manually  or  with  a computer)  it  is 
suggested  that  the  analyst  select  each  possible  set  of  two  extreme  points 

A 

in  y and  investigate  the  behavior  of  g(x  ; y)  as  a third  component  of 
y is  varied  and  X is  adjusted  to  maintain  optimality;  this  approach  is 
feasible  for  perhaps  three  or  four  observations.  For  more  observations, 
it  is  suggested  that  the  simple  linear  programs  ( 10)  and  ( 11)  be  solved 
directly  for  different  values  of  X 
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Unfortunately,  it  is  also  tedious  to  obtain  a confidence  region  for 

. *. 

g(X  ) alone.  The  values 


£>  = inf  g(x)  ( IT) 

\ € [_X  y X] 

g = supg(x)  (18) 

X e [X  , X] 

will  satisfy 

pr{_g  < g(x*)  <g)  > 1 - y (19) 

just  as  for  (9)  and  (12).  The  simplest  way  to  solve  ( 17)  and  (18)  is  by 
inspection  from  the  _g(x)  , 1(x)  plots. 

The  values  of  and  "g  do  not  in  general  occur  at  an  extreme  point 
of  y although  they  often  seem  to.  It  is  suggested,  therefore,  that  the 

^ /\  A 

U(y)  , g(x(y)  ; y))  values  be  plotted  for  each  feasible  y vector  with  three 
or  more  of  its  components  fixed  at  their  upper  or  lower  bounds.  This  will 
always  yield  the  values  for  \ and  X and  will  often  yield  the  values  for 
£ and  g . 
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J.  NUMERICAL  EXAMPLE 


I 
I 

^ In  this  section  we  present  a worked-out  example  as  an  aid  to  understand- 

- ing  the  analysis  of  the  preceding  section.  The  following  data  are  given: 


i 

ai 

A.1 

1 

A2 

1 

?i 

1 

-1 

1 

2 

.02 

2 

0 

k 

5 

.02 

3 

1 

1 

b 

.01 

Thus,  for  example,  [1,  2]  is  a 98/3  confidence  interval  for  g(x)  at 
X = -1  . The  probability  content  of  the  final  region  is  ( .98)(  ..98)(  .99)  = 

.95  approximately.  The  coefficients  are  given  by 


while  the  function  X(y)  is 


The  set  Y is  given  by  y = (y  | 1 < y^  < 2 , 4 < y0  < 5 , 1 < y . < 4) 

* 

The  bounds  on  X are  then  readily  found  by  observing  the  extreme  points 


of  the  set  y : 


X - 

10 


X = -± 


For  any  X e 


f.l  1 

L 10  ’ 2_ 


for  y = 


for  y = 


2 

b 

1 

any 

b 

b 


and  yeY  , we  have,  for  optimality, 


2X-  1 


y2  = 


b\ 


yl  + 


2X  - 1 

4x 


’ 


g(  X ; y)  = 


(1+X)2  (1-X)2 


4-X 


y3“  4x 


y1  U 4 o)  , 


g(0  ; y)  = y2  . 


The  resulting  g(x)  , "g(x)  functions  are  obtained,  in  this  case,  by 
examining  the  values  of  g(X  ; y)  when  two  components  of  y are  set  at  the 
various  extreme  points  and  the  third  component  is  allowed  to  vary.  Taking 
the  minimum  and  maximum  values  for  each  X yields: 


2X  +8X+4 

g(x)  = - ~ < X < \ (y  = 2 , y = 4) 

2X  + 1 ^ L * 


(so) 


12  - 


X1-  - 6x  + 1 


' 10  - } - 14  ^yl  2»  y3  “ ^ 


5(  >.)  = 


4x"  - lox  + 5 


1 - 2X 


4x"  + lox  + 5 


1 +2X 


1 U-XZ  0 


0 - x?  m 


(y2  = 5,  y3  = ^ 


(yx  = 1,  y2  = 5) 


3x  + lox  +3 


^<x<  \ (y1=i,  y =M 


The  segments  of  each  function  are  annotated  with  the  values  of  the  y vector 

components  which  remain  fixed.  The  functions  are  plotted  in  Figure  1.  Note 

the  characteristic  (and  unusual)  shape  of  the  joint  confidence  region  for 
* * . 

(X  , g(X  ))  . Figure  1 shows  the  upper  and  lower  confidence  bands  for  the 


response  g(X)  ; see  [3]  for  details. 
By  inspection  we  find  that 


_g  = 4 for  X = 0 

g = 5.225  for  X = 3/10 


for  numerical  example. 


4.  SOME  EXTENSIONS 


In  this  Section  we  take  up  a few  obvious  extensions  to  the  work  of 
Section  2.  These  are:  dealing  with  a quadratic  function  with  more  than 
three  simulation  outputs,  handling  non- concavity  and  selection  of  good  con- 
fidence interval  widths  at  each  base  point.  Extensions  to  the  numerical 
example  of  Section  3 are  also  presented. 

4. 1 Additional  Observations 

Suppose  that  a fourth  parameter  setting  a^  - l/2  is  chosen  for 
our  numerical  example  and  that,  with  probability  .98,  4 < g(  l/2)  < 5 

Including  this  additional  observation  decreases  the  probability  of  our  final 
results  to  about  ,93>  this  points  up  the  major  disadvantage  of  adding  new 
observations  to  an  existing  solution. 

The  various  optimization  problems  solved  in  Section  3 must  be  augmented 
by  adding  the  constraints 

h-~iyl+hy2+8y3~5  • 

When  this  is  done,  we  obtain 

X = - l/20  for  y = 

X = 1/2  for  y = 


2 

5 

4/3 

4 

2 

4 

4 

4 1/4 
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their  form  is  similar 


The  functions  _g(x)  and  g(x)  can  also  be  found; 
to  (20)  and  (21)  and  they  are  plotted  in  Figure  2.  The  results  for  the 
optimum  value  bounds  are 

£=41/2  for  X = l/4  , 

g = 5 l/6  for  X = l/4 

It  is  instructive  to  consider  the  reduction  in  size  of  the  various  con* 

* 

fidence  regions  when  the  fourth  observation  is  added.  The  intervals  for  X 
and  g(X  ) are  reduced  in  length  by  8.3$  and  8.2$,  respectively.  By  perform- 
ing a somewhat  tedious  integration,  one  may  find  that  the  joint  confidence 
region  on  (X  , g(X  ))  is  reduced  in  area  by  24$. 

4.2  Interval  Width  Selection 

Instead  of  performing  a simulation  with  an  entirely  new  parameter 
setting  as  in  the  previous  subsection,  one  may  choose  to  perform  more  repli- 
cations at  one  of  the  existing  settings  in  order  to  narrow  the  confidence 
1 2 

interval  [A^  , ] at  that  point.  Note  that  this  procedure  leaves  the 

probability  content  of  the  resulting  confidence  regions  unchanged.  Alter- 
natively, it  is  usually  possible  to  narrow  any  of  the  original  confidence 
intervals  if  one  is  willing  to  accept  the  resulting  decrease  in  the  proba- 
bility content  of  the  final  results.  We  thus  consider  the  selection  of  a 

1 2 

set  of  good  widths  for  the  [A^  , Af']  intervals. 

Because  of  the  complex  expressions  that  arise  if  one  allows  too  much 
generality,  it  is  difficult  to  give  robust  rules  for  choosing  the  interval 
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Figure  2.  Effect  of  adding  a fourth  base  point  to  example. 


widths.  Intuitively,  it  seems  best  to  select  all  with  about  the  same  width, 
or  else  with  about  the  same  probability  coverage.  In  order  to  consider  this 
matter  more  concretely,  we  shall  look  more  closely  at  the  example  of  Section  J. 


If  the  additional  simulation  replications  carried  out  at  X = ah  in 

Section  4.  1 had  been  run  instead  for  X = a^  , the  confidence  interval 
1 2 

[A,  , A^  ] would  be  based  on  twice  as  many  observations  and  would  thus  be 
about  70 . 7$  as  wide  as  the  original  interval  with  the  same  probability  coverage. 


The  calculations  of  Section  3 were  repeated,  therefore,  with  =1.5 

2 

A^  = 3.5  ; the  results  are  plotted  in  Figure  3 and  are  summarized  below; 


and 
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for 


y = 


2 

4 

1 2- 


= -5- 
14 


for 


y = 


1 

4 

3? 


£ = 4 


for 


X = 0 


g = 5 25//l76  for 


= 5 L 


X = '22 


The  size  of  [X  , X]  thus  decreases  by  31$  and  of  [_g  , g]  by  6.8$;  the  area 

*#-  -X- 

of  the  (X  , g(X  ))  region  decreases  by  41$.  It  may  thus  be  concluded  that 
in  this  instance  it  is  preferable  to  narrow  the  third  confidence  interval 
than  to  take  additional  observations  at  some  other  point. 


Figure  3-  Effect  of  narrowing  the  third  confidence  interval  in  the  example. 


To  pursue  the  matter  further,  we  investigate  how  the  widths  of  the 

— _ 12 
[X  , X]  and  [j*  , g]  confidence  intervals  change  as  the  [A^  , 1 interval 

changes.  We  assume  that 


A^1  = 2X/2  - 5 , = 2^2  + 5 . 

The  resulting  widths  are  plotted  in  Figure  4.  Note  that  when  5 = 3^2  in 


this  case,  the  vector  y = 14]  yields  a linear  function;  the  upper  bound 

\6  / 

X thus  approaches  + oo  as  5 increases. 

If  we  further  assume  that  the  simulation  response  at  X = a has  a 

3 

normal  distribution  (on  an  approximately  normal  distribution)  we  can  plot  the 
sizes  of  the  two  confidence  intervals  against  their  probability  coverage;  note 

that  since  the  a^  and  a^  intervals  do  not  change  in  this  case,  the  maximum 

2 

probability  is  (.98)  = .9604  . The  result  is  shown  in  Figure  5.  Note  the 

sharp  upward  turn  at  a point  corresponding  to  1 - 7 = .94  , i.e.  7^  ~ .02 

Both  Figures  4 and  5 lend  some  wieght  to  our  previous  observation  that 

* * 

the  "best"  confidence  regions  for  X and  g(  X ) result  when  the  absolute 
widths  of  the  simulation  output  confidence  intervals  are  about  the  same.  We 
must  qualify  this  statement  somewhat  when  the  parameter  settings  are  not  evenly 
spaced,  since  a greater  width  is  acceptable  for  parameter  settings  which  are 
relatively  distant  from  the  optimum. 

4.3  Non- Concave  Functions 

When  the  value  of  z (y)  found  from  (b)  is  not  negative,  the  result- 

3 

ing  g( X ) function  becomes  unbounded  as  X — ' i «>  ■ As  we  saw  in  Section  4.  _ , 
this  will  cause  either  X or  X or  both  to  become  unbounded  as  well. 
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is  normally  distributed. 


There  are  two  possible  types  of  response  to  the  situation  when  it 


turns  out  that  z^(y)  >0  . The  first  response  is  based  on  the  assumption 

that  g(x)  is  actually  a strictly  concave  function.  In  this  case  the 

analyst  attempts  to  adjust  the  simulation  results  so  as  to  obtain  concave 

functions.  The  other  type  of  response  merely  restricts  the  problem  ( l)  to 

the  closed  interval  [a^  , a ,]  where  & < a,  and  a . > a are 

O'n  + l 0—1  n + 1 — n 

user-specified  constants.  Of  course,  even  a convex  function  attains  a maxi- 
mum over  a closed  bounded  set.  In  the  rest  of  this  Section  we  discuss  these 
two  types  of  response  in  more  detail. 

The  value  of  z^( y)  will  be  negative  as  long  as 


(22) 


for  some  i<j<k  , n . The  adjustment  process  consists 

2 12 

of  modifying  one  triple  A.  , A.  , A,  if  none  of  them  satisfy  (22).  There 

i J x 

are  five  ways  to  do  this: 

( 1)  By  using  the  distribution  theory  that  was  used  to  obtain  the  original 

1 2 

[A^  , A^  ] intervals,  decrease  the  probability  coverage  (i.e., 
tighten  the  intervals)  until  (22)  is  satisfied.  It  may  also  be 
possible  (if  normal  theory  is  applicable)  to  use  asymmetric  rather 
than  symmetric  intervals  or  to  perform  some  other  adjustement.  It 
is  likely  that  only  a single  interval  will  need  to  be  altered. 

(2)  Carry  out  more  simulation  runs.  These  may  be  used  to  tighten  the 
intervals  already  obtained  or  to  investigate  the  response  at  an 
additional  base  point.  Based  on  the  results  of  this  section,  the 
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second  alternative  is  suggested  only  when  the  responses  at  the 
existing  points  indicate  that  the  optimum  lies  outside  the  current 
range. 

(3)  Inequality  (22)  may  be  expressed  in  terms  of  a constraint  on  y 

and  explicitly  introduced  into  the  various  optimization  problems 
of  Section  2 to  insure  that  z^(y)  < 0 . The  main  difficulty 

here  is  that  (22)  is  a strict  inequality;  for  this  reason,  a 
parameter  & should  be  chosen  and  the  following  constraint  used: 

a,-  aP  ap-  ai 

yP  > J — 2 y,  + — L y,  + <?  • (25) 

vai  vai 

A good  value  for  £ would  be  a few  percent  of  the  maximum 
response  obtained. 

(4)  The  work  of  Barlow,  et.  al  [1]  on  regression  under  order  restric- 

2 1 2 

tions  may  be  applied  to  adjust,  say,  , and  A^  so 

that  (22)  always  holds;  this  would  involve  estimates  of  the 
variability  (i.e.,  variance)  of  the  three  values.  This  approach 
allows  the  analyst  to  use  explicitly  the  additional  information 
that  g(x)  is  concave.  This  approach  is  somewhat  involved  and 
will  not  be  discussed  further  here. 

(5)  The  assumption  that  g(x)  is  concave  may  be  questioned  and  a 
statistical  test  performed  to  determine  its  likelihood.  This 
would  involve  information  on  the  distribution  of  the  simulation 
output  response,  however. 
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The  constraint  (23)  offers  a convenient  means  of  carrying  out  the  second 
approach,  i.e.  limiting  the  optimum  X to  the  interval  [a^  , a^ ^ . In 

this  case  we  solve  for  X and  X as  usual  (the  constraint  protecting  us 
from  catastrophic  unboundedness)  and  then  set 

X = maximum  ( a^  , _X) 

X = minimum  (a  , , X ) 
v n + 1’  ' 


As  long  as  $ is  quite  small  this  method  should  work  well. 
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